#ifndef _SYSTEM_HF_H_
#define _SYSTEM_HF_H_
#include<eigen3/Eigen/Eigen>//matrix library

#include"hfsolver_ho.h"

/*
a common interface for method based on but beyond HF method.
System_TzType and DataType can be changed according to needs.
*/

using Eigen::Matrix;
using Eigen::Dynamic;
using Eigen::SelfAdjointEigenSolver;

 
typedef double HFDataType;
/// system under HF basis
class System_Tz_HF:public System_Tz<SHF_Orbital<HFDataType>,HFDataType>
{
 public:
  typedef System_Tz_HO System_TzType;
  typedef SHFSolver_HO SHFSolverType;
  typedef typename System_TzType::DataType DataType;
  /* typedef System_Tz_WS System_TzType; */
  /* typedef SHFSolver_WS SHFSolverType; */
  /* typedef typename System_TzType::DataType DataType; */
  typedef Matrix<DataType,Dynamic,Dynamic> Mat;
  typedef Matrix<DataType,Dynamic,1> Vec;
  
  typedef SHF_Orbital<DataType> OrbitalType;
  const static int MatEledim;//< dimension of TwoBodyMatEle.
  typedef TwoBodyMatElement<DataType> TwoBodyMatEle;
  typedef Matrix<TwoBodyMatEle,Dynamic,Dynamic> Mat2B;

  System_Tz_HF(){}
 System_Tz_HF(System_TzType * Ptr):pSystem_Tz(Ptr),hfsolver(Ptr),hbar_omega(Ptr->hbar_omega){}
 System_Tz_HF(System_TzType * Ptr,int _A,int _Z):pSystem_Tz(Ptr),hfsolver(Ptr,_A,_Z),A(_A),Z(_Z),realA(_A),hbar_omega(Ptr->hbar_omega){}

  void setAZ(int _A,int _Z,int _realA=-1)
  {
    if(_realA==-1) realA=_A;
    else realA=_realA;
    A=_A;
    Z=_Z;
    hfsolver.setAZ(_A,_Z,realA);
  }

  void setupOrbitals();
  void setupTwoBodyMat();

  void setup(System_TzType * Ptr,int _A,int _Z)
  {
    pSystem_Tz=Ptr;
    hbar_omega=Ptr->hbar_omega;
    A=_A;
    Z=_Z;
    realA=_A;
    hfsolver.setup(Ptr,_A,_Z);
    setup();
  }
  void setup()
  {
    setupOrbitals();
    setupTwoBodyStates();
    setupTwoBodyMat();
  }
  ///interface for method beyond HF.
  DataType get1B(int bra,int ket) const;  
  DataType get2B(int a,int b,int c,int d,int J) const
  {
    TwoBodyMatEle MatEle;
    get2BmatAt(a,b,c,d,J,MatEle);
    //     get2BmatOnNeed(a,b,c,d,J,MatEle);
    if(MatEle.empty()) return 0;
    if(a==b) MatEle*=sqrt(2.0);
    if(c==d) MatEle*=sqrt(2.0);
    return MatEle[0]-MatEle[3]*pSystem_Tz->hbar_omega/realA;
  }

  ///interface for nocore shell model
  DataType get1B_p2(int bra,int ket) const;
  double get1B_p2_plus_r2(int bra,int ket) const;
  
  void get2BmatOnNeed(int a,int b,int c,int d,int J,TwoBodyMatEle & MatEle) const;

  void printOrbitals() const
  {
    hfsolver.printOrbitals();
  }

  
  ///maybe usefull operator
  ///only radial part
  DataType rL(int bra,int ket,int L) const;
  DataType rms() const
  {
    DataType val=0;
    for(int i=0;i<=FermiSurface;i++)
      {
	val+=rL(i,i,2);
      }
    return val/realA;
  }
  /*
    Q0=sum_i{r_i^2},S(E0)=2\hbar^2/m*A*<r^2>
  */
  ///be careful with the convention, Q0 is r^2, while QL proportion to r^L YL
  DataType Q(int p,int h,int J,int isospin) const;
  DataType Q1_v(int p,int h) const;
  DataType Q1_s(int p,int h,DataType Rms) const;
  
  int A,Z;//for HF. basis generating.
  int realA;//the real A in Hamiltonian which remove COM.
  int FermiSurface;
  DataType HFenergy;
  System_TzType * pSystem_Tz;
  SHFSolverType hfsolver;

  double hbar_omega;
 private:
  /// return <a|alpha>, a is the ath orbital of Orbitals in pSystem_Tz, alpha is alpha'th HF. orbital in Orbitals
  /* DataType overlap(int a,int alpha) */
  /* { */
  /*   if(Orbitals[alpha].GroupIndex != pSystem_Tz->GroupIndice[a]) */
  /*     return 0; */
  /*   return Orbitals[alpha].Coeff[pSystem_Tz->OrbIndice[a]]; */
  /* } */
};



#endif
